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SECTION  I 


INTRODUCTION 

The  ever  widening  application  of  computational  aerodynamics  to  complex 
engineering  problems  makes  the  need  for  efficient  numerical  algorithms 
more  apparent.  In  the  past  few  years,  several  numerical  algorithms  have 
been  developed'*",  but  additional  research  is  still  required  to  substantially 
improve  numerical  convergence  rate,  geometric  adaptability  and  reliability 
of  existing  methods.  Unfortunately,  no  attempt  has  been  made  to  establish 
a  commonly  accepted  criterion  for  the  evaluation  of  the  aforementioned 
requirements.  One  realizes  that  the  magnitude  of  this  undertaking  is  certainly 
beyond  the  scope  of  a  single  investigation.  However,  this  information  is 
vitally  important  to  the  application  of  computational  aerodynamics  in  the 
design  process.  Therefore,  the  evaluation  process  is  frequently  forced  to 
concentrate  on  the  clearly  defined  issue  of  relative  numerical  efficiency. 

In  the  present  effort,  a  comparative  study  was  focused  on  MacCormack’s 
explicit  scheme  and  its  implicit  analogue.  Specifically,  the  comparison  was 
carried  out  by  using  the  two  procedures  to  solve  a  practical  aerodynamic 
problem.  The  choice  of  solving  schemes  was  based  on  the  fact  that  the  explicit 
scheme  has  frequently  been  used  as  a  bench-mark"’ ’  ’  for  new  algorithm  develop¬ 
ment,  and  that  the  two  investigated  algorithms  are  easily  adapted  for  vector 
processing. 

The  numerical  simulation  of  flow  over  a  biconic  body  with  an  afterbody 
compression  flap  (Figs.  1,2)* was  obtained  in  this  study.  The  computation 
was  performed  for  flow  with  a  nominal  free  stream  Mach  number  of  7.97  and 
a  characteristic  Reynolds  number  of  9.23  million.  This  particular  geometric 
configuration  generates  a  wide  variety  of  aerodynamic  characteristics  not 

*  Figures  begin  on  page  24. 
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usually  encountered  in  a  single  configuration.  In  particular,  the  flow 
field  contains  an  enveloping  bow  shock,  rapid  expansions  at  the  juncture 
of  the  biconic  forebody  and  the  sliced  conical  afterbody  and  finally  an 
inviscid-viscous  interaction  around  the  afterbody  compression  flap.  The 
solving  procedure  also  offers  the  opportunity  to  evaluate  the  idea  of 
segmented  and  overlapped  computational  domains^ *^.  Basically,  the  problem 
was  separated  into  two  interdependent  computational  domains.  The  forebody 
was  simulated  by  a  simple  axisymmetric  numerical  formulation.  The 
comparative  study  of  the  explicit  scheme  and  its  implicit  analogue  was 
also  limited  to  the  axisymmetric  forebody  region.  The  more  complex  flow 
field  around  the  afterbody  was  then  analyzed  using  the  full  three-dimensional 
equations.  The  solution  of  the  complete  configuration  is  presented  in 
composite  form.  The  specific  comparison  with  experimental  data^  in  terms 
of  surface  static  pressure  and  heat  transfer  distributions  as  well  as  pitot 
pressure,  Mach  number  and  tangential  velocity  profiles  will  be  given  to 
verify  the  accuracy  of  the  numerical  simulation.  The  entire  flow  structure 
around  the  biconic  configuration  with  the  afterbody  compression  flap  in 
terms  of  density  contours,  cross-flow  velocity  distribution  and  surface 
shear  stress  mapping  will  also  be  depicted.  Finally,  a  detailed  presentation 
on  the  comparative  study  of  the  explicit  scheme  and  its  implicit  analogue 


will  be  delineated. 


SECTION  II 


ANALYSIS 


The  time  dependent,  three-dimensional  Navier-Stok.es  equations  in 
mass-averaged  variables  and  in  the  chain-rule  conservation  law  form  can 
be  given  as^ 


The  dependent  variables  are  U(p,  pu,  pv,  pw,  pe) .  The  flux  vectors,  F,  G 
and  H  are  simply  the  Cartesian  components  of  the  continuity,  momentum  and 
energy  equations.  The  axisymmetric  formulation  is  but  a  subset  of  the 
governing  equations. 
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(4) 


where,  the  component  of  shear  stress  is  defined  by 
Ti^  =  (p+e)  (Def  u).^  -  [-|(p+e)V*u+p]  6^ 


(5) 


The  closure  of  this  system  of  equations  is  achieved  by  introducing  the 

9  8 

Baldwin-Lomax  turbulence  model  with  a  minor  modification  and  by  assigning 
a  turbulent  Prandtl  number  of  0.9.  Specifically,  the  two-layer  eddy  viscosity 
model  is  given  by 
Inner  region: 

e  =  p (0.4LD)2  | ml  (6) 

In  the  present  formulation,  the  vorticity  of  the  flow  w,  the  Van  Driest 


damping  factor  D,  and  the  scaling  length  L  are  given  as  follows 
1  „ 

u  =  v  x  u 

P 1  I  I  T  /  9 

D  =  1-exp  [ - — - )  L/26] 


(7) 

(8) 


l  *  i(x-xb)2  +  (y-yb)2  +  (z-Zjj)2]1^2 


Outer  region: 

e  =  0.0336  pF  / [1+5. 5 (0. 3L/L  )6] 
waxe  max 


(9) 


(10) 


The  wake  function,  Fwake  is  the  minimum  value  of  the  two  following 
expressions  at  any  point  in  space. 

(11) 


F  =  L  *  F 
wake  max  max 

F  ,  =  0.25  •  L  •  (u2+v2+w2)  /F 

wake  max  max  max 


(12) 


In  the  present  analysis,  is  the  value  of  the  length  scale  where 

F=LD|w|  reaches  it’s  maximum  value  within  the  turbulent  shear  layer.  The 
system  of  equations  is  formally  closed  with  Sutherland’s  equation  for  molecular 
viscosity  and  the  equation  of  state  for  a  perfect  gas. 
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Since  a  composite  solution  obtained  by  overlapping  computational 
domains  for  the  investigated  problem  was  performed,  different  initial 
and  boundary  conditions  were  imposed  for  each  subproblem.  However,  since 
these  conditions  are  quite  similar,  they  will  be  discussed  together.  The 
initial  condition  for  the  forebody  solution  consisted  of  freestream  values 
for  all  field  points  excluding  the  surface  nodes. 

U  (0,5, n,C)  =  0^  (13) 

The  initial  condition  of  the  three-dimensional  afterbody  solution  is  described 
in  part  by  the  forebody  calculation.  The  overlapping  station  is  located  at 
a  distance  of  45.5  radii  downstream  of  the  blunt  nose.  The  initial  value 
assigned  to  the  overlapping  surface  also  serves  as  the  upstream  boundary 
condition  for  the  afterbody  calculation.  This  is  the  only  link  between  the 
forebody  and  the  afterbody  solutions. 

On  the  solid  body  contour,  the  no-slip  condition  for  velocity  components, 
the  isothermal  condition  for  temperature  to  duplicate  experimental  conditions, 
and  the  ortho-isobar ic  condition  for  pressure  are  imposed.  The  value  of 
surface  density  is  deduced  from  the  equation  of  state. 

u,  v,  w  =  0  (14) 

Tb  =  311. 3°K  (15) 

n • wp  =  0  (16) 

where  the  outward  normal  of  the  body  oriented  coordinate  is  simplv, 
n=Vn/ | | Vn | | . 

Since  the  detached  bow  shock  wave  isolates  the  interacting  flow  domain, 
the  far  field  condition  requires  that  the  flow  remain  unperturbed. 

u  (t,e,i,o  =  (i7) 

For  the  afterbody  computation,  the  present  analysis  takes  advantage  of  the 
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property  of  symmetry  with  respect  to  the  y  axis  (Fig.  3).  Only  half  of 
the  cross-flow  plane  is  evaluated.  A  reflection  condition  is  applied 
which  insists  that  the  c,  component  of  velocity  be  equal  to  zero  on  the 
plane  of  symmetry. 

an 

—  *  0,  w  =  0  (18) 

Finally,  the  usual  no-change  condition  is  imposed  at  the  downstream  boundary 
for  each  computational  domain.  This  system  of  boundary  conditions  is  known 
to  be  well-posed  and  stable  for  the  supersonic  problem. 


SECTION  III 


COORDINATE  SYSTEM  AND  GRID  GENERATION 


The  numerical  generation  of  body-oriented  coordinates  is  straight- 
jnrward.  The  adopted  methods  is  the  homotopy  scheme^.  Basically,  it 
is  an  interpolation  procedure  between  grid  control  surfaces.  Once  the 
inner  and  the  outer  surface  are  given,  any  number  of  mesh  surface  points 
is  controlled  by  an  exponential  stretching  function.  For  turbulent  flow 
this  stretching  is  nearly  optimum  in  describing  the  logarithmic  velocity 
profile.  The  mesh  formulation  is  given  by 


.  ,,  ekri-l, 

y  ’  yo(_k  r*  +  vi (l-  - — > 


e  -1 


k  , 
e  -1 


.ekn-l1  ,  ..  ekn-l* 

Z  =  Z  (“5 - )  +  Z.  (1 - : - ) 

ok.  i  k 


(19) 


(20) 


e  -1  e  -1 

where  0<nj[l .  The  exponent  of  the  stretching  function  is  uniquely  determined 
by  the  minimum  distanced  between  the  body  and  the  immediately  adjacent 
coordinate  surface,  and  by  the  number  of  points  used.  Additional  mesh  clus¬ 
tering  around  the  bow  shock  wave  is  generated  by  increasing  the  total  number 
of  nodes  in  n  direction  by  four  in  this  region  and  subdividing  the  mesh  pre¬ 
viously  generated  using  the  logarithmic  coordinate  surfaces  locally  near  the 
bow  shock  (Figure  2  and  3).  For  the  baseline  case  of  computation,  the  calculated 
value  of  the  exponent  k  is  4.820.  The  finest  outward  normal  distanced  from 
the  body  surface  is  assigned  a  value  of  0.00305  cm  in  the  blunt  nose  region. 
However,  the  grid  spacing  increases  four-fold  in  the  aft  portion  of  the  biconic 
forebody,  where  the  law  of  the  wall  variable  y+  assumes  a  value  of  5.901. 

For  the  forebody  calculation,  the  outer  coordinate  line,  0*1,  is  described 
by  a  curve  consisting  of  a  circle  and  a  parabola  matched  to  closely  approximate 
the  bow  shock  wave.  Care  is  also  exercised  to  ensure  that  the  generated 
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coordinates  are  nearly  orthogonal.  In  order  to  achieve  a  systematic 
comparative  study,  various  mesh  distributions  are  generated  by  reducing 
the  An  value  (inversely  proportional  to  the  total  number  of  nodes  chosen) 
and  retaining  the  exponent  k  at  the  fixed  value  of  4.820.  In  this  manner, 
the  coordinate  systems  are  affine  and  the  computational  domains  are  unaltered. 
The  streamwise  numerical  resolution  is  also  evaluated  by  using  two  different 
streamwise  mesh  systems  which  comprise  (62x40)  nodes  and  (90x40)  nodes, 
respectively. 

The  three-dimensional  afterbody  calculation,  the  coordinate  system 
is  constructed  by  a  series  of  consecutive  axial  cross  sections  evenly  spaced 
to  achieve  optimal  numerical  resolution.  The  first  thirteen  streamwise 
step-spacings  upstream  of  the  compression  flap  have  the  value  of  0.859  cm, 
which  is  less  than  half  the  boundary-layer  thickness  at  the  overlapping 
plane.  A  finer  streamwise  step  size  (0.4276  cm)  is  used  for  the  rest  of 
the  streamwise  stations.  The  body  geometry  variation  is  relatively  mild 
with  two  sliced  surfaces  on  the  top  and  bottom  of  the  conic  body  and  a  20° 
compression  flap.  Therefore,  only  38  circumferential  points  were  used  to 
define  the  cross-sections.  All  the  mesh  points  have  an  even  angular  dis¬ 
placement  of  5°.  The  first  two  and  the  last  two  arrays  of  mesh  points  were 
astriding  the  plane  of  symmetry.  A  total  of  (33,  40,  38)  mesh  points  were 
used  for  the  afterbody  calculation.  In  Figure  3,  four  typical  cross-sectional 
grid  point  distributions  are  presented. 


SECTION  IV 


Since  this  stability  condition  for  time  step  size  does  not  contain  viscous 
terms,  the  maximum  allowable  time  increment  for  the  actual  computation  must 
be  determined  numerically.  For  the  present  analysis,  the  highest  and  most 
consistent  CFL  number  is  0.8  and  is  used  throughout  as  the,  criterion  for 
the  numerical  efficiency  study. 

The  detailed  description  of  MacCormack's  implicit  algorithm  is  documented 
in  reference  3  and  therefore  will  not  be  elaborated  on  here.  The  implicit 
algorithm  can  be  directly  implemented  into  the  existing  code,  with  the  basic 


idea  being  to  convey  the  locally  determined  solution  evolution  globally 
in  a  stable  manner.  The  implicit  propagation  of  local  solution  changes 
in  time  are  given  by  the  following  equation. 


.  .  3A*  ,  3Bv  3Un+1  3Un 

(I+At  3n  At  3c  ^  3t  *  3t 


(25) 


where  A  and  B  are  the  Jacobians  of  F  and  G  with  respect  to  the  dependent 
variable  U.  The  notation  in  the  (A*,  B*)  above  expression  simply  indicates 


3Un+1 

that  the  spatial  derivative  also  operates  on  the  time  derivative  — ^ —  . 

The  numerical  procedure  is  first  performed  for  the  explicit  predictor  sweep, 
then  followed  by  the  implicit  sweep  for  the  nodes  where  the  time-step 
exceeds  the  CFL  limit  locally.  This  procedure  is  repeated  for  the  corrector 
sweep.  The  implicit  differencing  operator  is  parallel  to  that  of  the  explicit 
differencing  operator.  The  alternative  one-sided  difference  approximates 
the  first  derivative  and  the  second  derivative  are  centrally  differenced. 


AJa|  a  Ib| 

(I-At  )(I-At  —77 — )  6Un+1  =  AUn 


AC 


(26) 


A_ 1 A I  A  |B| 

(I+At  -77— )(I+At  -7 — )  6Un  1 


An 


AC 


AU 


(27) 


where 


n+1 


At 


3Un+1 

3t 


AUn 


At 


3Un 

3t 


The  matrices  [ A |  and  | B |  are  the  matrices  with  positive  eigen  values  and 
are  related  to  the  Jacobians  A  and  B.  Specifically,  the  matrices  A  and  B 
can  be  expressed  by  it’s  diagonalized  matrices  and  the  similar  transformation 
matrices  in  terms  of  the  general  coordinates 
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The  flux  change  in  either  the  n  or  the  ?  directions  are  easily  given  as 
s'1  3U,  and  S^1  3U,  respectively. 

Since  the  procedure  is  implicit,  the  treatment  of  boundary  conditions 
is  critical  to  further  improve  the  rate  of  numerical  convergence.  At  this 
writing,  the  transmitted  flux  change  at  the  £  *  0  boundary  is  swept 
immediately  back  into  the  flow  field,  leaving  no  corrections  to  be  performed 
later.  However,  the  present  procedure  stored  the  transmitted  flux  received 
at  the  inner  boundary  (n  =  0)  during  the  predictor  implicit  sweep,  and  then 
later  sweeps  it  back  out  during  the  corrector  step.  Although  flux  changes  are 
rigorously  conserved,  there  is  a  time  lag.  The  implementation  of  implicit 
boundary  condition  is  a  vital  area  for  future  research  and  continuous  efforts 
are  currently  still  being  undertaken. 


For  the  investigated  case,  the  static  pressure  jumped  seventy-four  fold 
across  the  bow  shock  wave  enveloping  the  blunt  nose  region.  Numerical 
damping  was  required  to  surpress  numerical  oscillations,  therefore,  fourth- 
order  pressure  damping  terms'1’1  are  used  in  the  present  analysis.  The 
magnitude  of  the  artificial  viscosity-like  terms  are  restricted  to  a  value 
not  to  exceed  0.5. 

At  An3  UuJ+(n2+n2+n2)1/2C]  ^  |-Mr|  <0.5 


A  transient  damping  term  representing  a  measure  of  the  change  in 

entropy  primarily  induced  by  viscous  dissipation  is  also  used  to  maintain 
3 

stability  . 

|6p/c2  -  6p|/(|t)  (izi)  p 

This  additional  precautionary  measure  vanishes  when  the  asymptotic  steady 
state  solution  is  reached. 

The  aforementioned  implicit  numerical  procedure  results  in  matrix 
equations  to  be  solved  that  are  either  upper  or  lower  block  bidiagonal. 

Since  no  recurrent  relationship  is  required  to  invert  the  system  of  equations, 
the  vectorization  of  the  implicit  procedure  is  easily  achieved.  The 
implicit  calculation  needs  only  twenty-three  percent  more  CPU  time  to  process 
the  data  than  the  explicit  scheme.  The  additional  computing  resources 
required  for  the  implicit  code  are  nearly  identical  for  either  scalar  or 
vector  processors.  Finally,  the  numerical  convergence  criterion  for  the 
present  analysis  is  defined  as  a  change  in  the  Stanton  number  of  less  than 
one  percent  per  one  half  of  a  characteristic  time.  The  time  scale  is  the 
period  for  a  fluid  particle  to  travel  from  the  upstream  to  downstream  boundaries 
of  the  computational  domain  at  the  freestream  velocity. 


SECTION  V 


DISCUSSION  OF  RESULTS 

The  present  results  will  be  discussed  in  two  parts.  In  the  first 
topic,  the  comparative  study  of  the  explicit  method  and  it's  implicit 
analogue  applied  to  simulate  the  biconic  forebody  is  presented.  A  brief 
discussion  is  also  included  on  the  relative  numerical  efficiency  between 
the  full  and  the  parabolized  Navier-Stokes  solutions  for  the  investigated 
configuration.  The  rest  of  the  presentation  concentrates  on  the  delineation 
of  the  detailed  flow  structure  around  the  afterbody  with  a  20°  compression 
flap  and  a  sliced  conic  shape.  All  of  the  comparisons  and  validation  are 
accompanied  by  experimental  data  collected  under  identical  conditions. 

The  comparative  study  of  MacCormack's  explicit  scheme  and  it's  implicit 
counterpart  is  focused  on  the  axisymmetric  forebody  from  the  following  ration¬ 
al.  First,  the  major  attraction  of  using  an  implicit  scheme  to  simulate  a 
given  problem  is  that  of  the  favorable  stability  characteristic  of  the  implicit 
scheme  leading  to  a  decrease  in  the  required  computing  resources.  In  principle, 
the  implicit  scheme  always  requires  more  arithmetic  operations  than  the  explicit 
method  to  advance  a  given  time  step.  Thus  the  favorable  allowable  time  step 
size  must  overwhelm  the  conditionally  stable  explicit  scheme  to  gain  an  overall 
advantage  in  data  processing  time.  For  most  aerodynamics  applications,  the 
mesh  refinement  for  the  viscous  dominated  region  usually  is  not  required  along 
more  than  two  coordinate  directions*’5’8’9.  Therefore,  the  implicit  procedure 
may  not  need  to  be  implemented  for  all  the  coordinates  in  order  to  achieve  the 
most  efficient  results.  Second,  the  forebody  numerical  simulation  offers  a 
severe  challenge  for  any  numerical  algorithm  in  simulating  the  strong  bow 
shock  wave  around  the  rapidly  changing  curvature  of  the  blunt  nose  and  the 


resulting  rapid  flow  expansion  that  follows.  Finally,  each  transformed 
coordinate  must  rotate  90°  to  accomodate  the  definition  of  the  body. 

In  the  following  table,  the  data  processing  rate  and  the  total  CPU 

time  used  to  meet  the  convergence  criterion  mentioned  earlier  are  given.  It 

is  observed  that  the  ratio  of  data  processing  rates  between  the  explicit 

scheme  and  the  implicit  scheme  remains  nearly  identical,  regardless  of  whether 

the  codes  were  performed  on  scalar  or  vector  processors.  This  is  a  clear 

indication  that  MacCormack's  implicit  scheme  is  more  readily  vectorized  than 

other  implicit  methods  which  require  tridiagonal  block  matrix  inversion. 

This  advantage  in  application  is  critical  for  large  scale  data  processing 

at  present.  The  allowable  time  step  size  for  the  implicit  algorithm  is 

recorded  for  a  stringent  application  standard;  the  accepted  CFL  number  must 

be  sustained  until  the  evolving  solution  meets  the  convergence  criterion. 

No  prolonged  relaxation  period  was  allocated  to  diminish  the  larger  temporal 

truncation  errors.  Most  important  of  all  only  the  practical  fine  grid 

spacing  adjacent  to  the  surface  sufficient  to  resolve  the  heat  transfer  was 

used.  In  the  eddy  viscosity  formulation,  the  laminar  sublayer  thickness 

is  normalized  by  the  law  of  the  wall  variable  y+  with  a  value  of  26.  Based 

on  early  investigations  of  strong  invisc id-viscous  interactions,  a  y+  value 

less  than  10  is  usually  adequate  to  resolve  the  shear  stress  and  heat  transfer 
12  13 

information  ’  .  Further  refinement  of  the  mesh  adjacent  to  the  surface 

would  yield  a  higher  CFL  value,  but  may  not  be  cost  effective.  Nevertheless, 
for  the  case  of  y  =  5.901,  the  implicit  procedure  can  consistently  process 
data  at  a  time  step  value  A. 125  times  higher  than  the  explicit  scheme.  The 
net  total  computing  resources  savings  amounts  to  a  factor  of  3.3  in  comparison 
with  the  explicit  procedure.  As  the  surface  mesh  spacing  reduces  to  correspond 


to  a  y+  value  of  3,  the  implicit  scheme  performs  at  a  consistent  CFL 
value  of  16.1.  The  computational  resource  savings  increases  to  an 
impressive  factor  of  10.09.  The  difference  in  predicting  the  Stanton 
number  distribution  between  the  fine  (y+  =  3)  and  the  course  (y+  -  5.9) 
mesh  solution  is  less  than  one  percent.  For  our  purpose,  the  numerical 


results  are  considered  identical. 


Algorithm 

Mesh  Spacing 

n 

Computer 

DPR 

_ 

Total  CPU 

At 

AtCFL 

Explicit 

(62,40) 

5.9 

CYBER  175 

6.1x10  A 

15123.0 

1.0 

Implicit 

(62,40) 

5.9 

CYBER  175 

7.5x10"* 

4582.7 

4.125 

Explicit 

(90,40) 

5.9 

CRAY  IS 

1. 78xlO~5 

641.9 

1.0 

Implicit 

(90,40) 

5.9 

CRAY  IS 

2.21xl0“5 

194.5 

4.110 

Explicit 

(62,40) 

5.9 

CRAY  IS 

1.72xl0"5 

426.00 

1.0 

Implicit 

(62,40) 

5.9 

CRAY  IS 

2. 13xl0~5 

129.10 

4.110 

Explicit 

(62,41) 

3.0 

CRAY  IS 

1. 72xl0~5 

1864.2 

1.0 

Implicit 

(62,41) 

3.0 

CRAY  IS 

2.43xl0"5 

142.0 

16.1 

The  afterbody  calculation  was  accomplished  using  a  (33,40,48)  mesh  system  to 
describe  the  flow  field  around  the  20°  compression  flap.  The  initial  condi¬ 
tion  for  the  afterbody  simulation  was  prescribed  by  the  value  of  the  over¬ 
lapping  plane,  which  is  a  much  better  guess  than  the  forebody  initial  condi- 
tion.  At  a  data  processing  rate  of  6.1x10  seconds,  the  numerical  solution 
converged  after  2,374.8  seconds  CPU  time.  On  the  other  hand,  the  numerical 
solutions  generated  by  various  parabolized  Navier-Stokes  codes'*^’16  of  a 
similar  configuration  (10°  flap)  required  a  range  of  291  to  537  seconds  of 
computer  time  on  a  CARY  computer16.  The  numerical  efficiency  of  parabolized 
Navier-Stokes  procedures  is  impressive.  However,  the  cited  CPU  time  excludes 
the  Navier-Stokes  calculations  of  the  blunt  nose  which  these  procedures  need 
as  an  initial  condition.  The  difference  in  total  computing  resources  for  the 
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entire  numerical  simulation  between  the  full  Navier-Stokes  and  the  para¬ 
bolized  Navier-Stokes  equations  is  less  than  an  order  of  magnitude. 

In  Figure  4,  the  comparison  of  the  surface  pressure  distributions 
on  the  forebody  reveals  excellent  agreement.  The  solutions  by  the  explicit 
and  implicit  methods  are  indistinguishable.  The  maximum  discrepancy  between 
experimental  data  and  calculations  is  limited  to  less  than  four  percent. 

The  Stanton  number  distributions  and  the  accompanying  data  is  given  in 
Figure  5.  The  over-all  agreement  between  data  and  the  calculations  is  very 
good.  The  maximum  discrepancy  is  near  the  high  expansion  and  recompression 
region  where  laminar  turbulent  transition  might  occur.  Since  a  tripping 
device  was  used  for  the  experimental  measurement  to  ensure  a  fully  developed 
turbulent  flow,  the  numerical  simulation  would  not  be  expected  to  duplicate 
the  local  phenomenon.  However,  downstream  of  the  junction  of  the  biconic 
configuration,  agreement  between  the  experimental  data  and  the  calculation 
is  excellent.  The  prediction  by  the  Navier-Stokes  equations  indicates  a 
comparable  result  to  that  of  the  parabolized  Navier-Stokes  equations1^.  In 
this  comparison,  the  implicit  solution  reveals  merely  a  1.6  percent  higher 
value  in  Stanton  number  evaluation  than  the  explicit  result.  This  is 
also  the  maximum  disparity  between  the  implicit  and  the  explicit  solutions 
investigated . 

The  comparisons  of  pitot  pressure  and  tangential  velocity  profiles 
are  given  in  Figures  6  and  7,  respectively.  The  detailed  flow  field  survey 
was  performed  at  41.5  radii  downstream  of  the  blunt  nose  in  the  meridian  plane. 
The  agreement  between  experimental  data  and  the  calculations  is  very  good. 
Again,  differences  between  the  solutions  achieved  by  the  explicit  method 
and  the  implicit  scheme  are  indistinguishable.  The  major  deviation  between 
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the  experimental  data  and  the  numerical  solutions  is  confined  within 
the  entropy  layer  induced  by  the  bow  shock  curvature  in  the  stagnation 
region,  and  in  the  definition  of  the  shock  wave.  This  is  a  fundamental 
short  coming  of  a  shock  capturing  scheme. 

In  order  to  better  organize  the  discussion  of  the  three-dimensional 
afterbody  calculation,  the  density  contour  of  the  entire  biconic  configura¬ 
tion  with  a  20°  compression  flap  in  the  meridian  plane  is  depicted  in 
Figure  8.  The  enveloping  bow  shock  wave,  the  rapid  expansions  downstream 
of  the  stagnation  region  and  the  conic  junction  and  finally  the  coalesced 
compression  shock  waves  over  the  flap  are  clearly  indicated.  In  the  lower 
half  of  the  meridian  plane  at  the  end  of  the  conical  section,  a  continuous 
expansion  is  also  observed  over  the  bottom  sliced  surface.  To  give  a  proper 
perspective  of  the  three-dimensional  flow  field  structure,  four  cross- 
sectional  density  contours  are  given  in  Figure  9.  The  systematic  development 
of  the  expansion  zone  over  the  slice  and  compression  flap  is  obvious.  The 
most  interesting  feature  is  that  the  coalescing  shock  wave  system  over  the 
flap  is  isolated  by  the  expansion  domain  above  the  slice.  Therefore,  the 
distortion  of  the  envelopiig  bow  shock  is  rather  limited. 

The  surface  pressure  distributions  on  the  body  surface  in  the  meridian 
plane  are  presented  in  Figure  10.  The  agreement  on  both  the  upper  and  the 
lower  surfaces  are  excellent.  A  good  agreement  also  is  reached  for  the 
circumferential  Stanton  number  distributions  between  data  and  calculation 
in  two  streamwise  locations,  when  experimental  measurements  are  available 
(Figure  11).  The  angular  displacement,  defined  to  be  zero  in  the  upper 
meridian  plane,  increases  in  the  clockwise  direction  until  the  lower 
meridian  plane  is  reached.  The  streamwise  Stanton  number  distribution  is 
also  calculated  and  the  agreement  with  data  is  very  good. 


Four  pitot  pressure  surveys,  one  upstream  of  the  location  of  the 
compression  flap  (x/rN  =  49.01),  one  roughly  one  boundary-layer  thickness 
downstream  of  the  flap  junction  and  two  locations  near  the  end  of  the 
compression  flap  are  illustrated  in  Figure  12.  The  last  survey  data  set 
was  collected  in  the  model  centerline  and  near  the  spanwise  edge  of  the 
compression  flap.  Very  good  agreement  between  data  and  the  calculation  is 
observed.  The  maximum  deivation  is  estimated  to  be  confined  within  a  few 
percent  and  is  mostly  limited  to  the  domain  of  the  entropy  layer  and  the 
definition  of  the  bow  shock.  The  peak  values  of  the  pitot  pressure  around 
the  coalescing  shock  wave  are  predicted  with  increasing  accuracy  in  the 
fine  mesh  distribution  near  the  body  surface. 

The  surface  shear  over  the  entire  afterbody  is  presented  in  Figure  13. 

The  0=0  ordinate  indicates  the  upper  meridian  plane  and  the  0  =  180° 
represents  the  lower.  This  picture  reveals  the  flow  structure  in  remarkable 
clarity  over  the  two  sliced  surfaces  where  the  flowstream  expands  rapidly 
toward  the  meridian  plane.  However,  the  flow  over  the  conic  section 
(45°<0_<135O)  possesses  no  swirling  motion.  On  top  of  the  upper  sliced 
surface,  which  is  followed  by  the  20°  compression  flap,  the  expansion  outward 
to  the  conic  section  of  the  body  and  the  expansion  inward  to  the  meridian 
plane  stands  out.  Three-dimensional  phenomena  are  clearly  evident.  Reversed 
flow  is  indicated  at  the  junction  of  the  sliced  surface  and  the  compression 
flap.  In  an  enlarged  streamwise  velocity  distribution  plot  in  the  meridian 
plane  (Figure  14),  the  reversed  flow  is  embedded  in  an  extremely  small  domain. 

In  the  last  figure.  Figure  15,  the  cross-flow  velocity  distribution  due 
to  the  configuration  asymmetry  is  presented.  The  cross-flow  velocity  com¬ 
ponents  are  given  in  radial  and  tangential  components.  This  picture  not 


SECTION  VI 


CONCLUSIONS 

The  numerical  simulation  of  flow  over  a  biconic  body  with  a  compression 
flap  is  successfully  accomplished  at  a  Mach  number  of  7.97,  a  Reynolds 
number  of  9.228  million  and  at  a  zero  degree  angle  of  attack.  The  numerical 
results  duplicate  all  the  experimental  measurements  in  pressure,  pitot 
pressure  and  heat  transfer  within  a  few  percent.  The  numerical  results 
simulate  the  intricate  flow  field  structure  around  a  three-dimensional 
configuration. 

The  application  of  the  MacCormack's  implicit  algorithm  indicates  the 
procedure  is  easily  vectorized  for  large-scale  data  processing,  and  a 
significant  gain  in  computer  resources  is  realized.  However,  continuous 
research  efforts  to  achieve  a  greater  numerical  efficiency  are  still  needed. 
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Figure  3  Mesh  System  of  the  Biconic 
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Figure  5  Comparison  of  Heat  Transfer 
on  Forebody 
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Figure  6  Comparison  of  Pitot  Pressure 
Profile,  X/Rn  =  41.45 


TANGENTIAL  VELOCITY  PROFILE 


Figure  7  Comparison  of  Tangential 

Velocity  Profiles  (X/Rn  =  41.43) 


Figure  9  Density  Contours  in  Cross- 
Sectional  Planes 


SURFACE  PRESSURE  DISTRIUTION 


Figure  10  Surface  Pressure  Distributions 
in  Meridian  Plane 


Figure  11  Comparison  of  Circumferential 
Stanton  Number  Distributions 
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Figure  12  Comparisons  of  Pitot  Pressure 
Profiles  (Continued) 


Comparisons  of  Pitot  Pressure 
Profiles  (Concluded) 


VELOCITY  PROFILE  IN  MERIDIAN  PLANE 


Figure  14  Velocity  Field  Over  the 
Compression  Flap 
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Figure  15  Cross-Flow  ' 
tion,  (X/Rn 
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